library(scuttle)
library(BiocParallel)
library(SingleCellExperiment)
source("code/_utils.R")
bp <- MulticoreParam(3)
set.seed(251103)jst
preamble
loading
pbs <- readRDS("outs/qbs.rds")
sce <- readRDS("outs/fil.rds")
ist <- readRDS("outs/lv1.rds")
lapply(pbs, colnames)$mye
[1] "DC" "mono.c" "mono.nc" "pDC"
$nkt
[1] "NK/ILC" "Tc" "Tcm" "Tcn" "Tex" "Tfh" "Thn" "Tp"
[9] "Treg"
sub <- sort(unique(ist$clust))
(names(sub) <- sub)[1] "mye" "nkt" "str" "tum"
selection
names(ks) <- ks <- unique(colnames(mx <- ist$profiles))
mx <- sapply(ks, \(k) rowMeans(mx[, k, drop=FALSE]))
sel <- lapply(ks[ks != "tum"], \(.) {
my <- mx[, setdiff(ks, .)]
fc <- rowMins(apply(my, 2, \(x) mx[, .]/x))
fc <- fc[!is.na(fc) & is.finite(fc) & fc > 1]
gs <- names(tail(sort(fc), 2e3))
})
sapply(sel, length)str mye nkt
841 729 360
clustering
jst <- bplapply(sub[sub != "tum"], \(.) {
ns <- c(1e4, 2e4, 1e5)/4
nk <- seq(2, ifelse(. == "mye", 4, 8))
if (is.null(gs <- sel[[.]])) gs <- TRUE
se <- sce[, names(which(ist$clust == .))]
ist <- .ist(se, nk=nk, ns=ns, gs=gs, pbs=pbs[[.]])
ks <- intersect(letters, unique(ist$clust))
ks <- c(ks, sort(colnames(pbs[[.]])))
ist$clust <- factor(ist$clust, ks)
ist
}, BPPARAM=bp)kst <- ist
idx <- ist$clust == "tum"
kst$prob <- kst$prob[idx]
kst$logliks <- kst$logliks[idx, ]
kst$clust <- factor(kst$clust[idx])
jst$tum <- kstvisuals
Code
lapply(setdiff(sub, "tum"), \(.) {
kid <- droplevels(kid <- jst[[.]]$clust)
sce$ist <- kid[match(colnames(sce), names(kid))]
.plt_fq(sce, "sid", "ist", id=., h=TRUE)
}) |> wrap_plots(nrow=1)Code
ps <- lapply(setdiff(sub, "tum"), \(.) {
kid <- droplevels(kid <- jst[[.]]$clust)
sce$ist <- kid[match(colnames(sce), names(kid))]
.plt_fq(sce, "ist", "sid", id=.)
})
ws <- sapply(ps, \(p) nlevels(p$data$Var1))
wrap_plots(ps, guides="collect", widths=ws)[[1]] NULL
[[2]] NULL
[[3]] NULL
labeling
lab <- list(
str=list(
epi="a",
FRCcts="b",
BEC="c",
LEC="d",
FRCtcz="e",
FRCpv="f",
fib="g",
mCAF="h"),
mye=list(
tum="a",
mac.pi="b",
AML="d",
DC="DC",
pDC="pDC",
TAM="c",
mono.c="mono.c",
mono.nc="mono.nc"),
nkt=list(
Tha="a",
`NK/ILC`="b",
Tp="Tp",
Tex=c("Tex", "c"),
Tfh="Tfh",
Treg="Treg",
Tcn="Tcn",
Tcm="Tcm",
Thn=c("Thn", "d"),
NK="NK/ILC"))
kst <- lapply(names(lab), \(.) .jst(jst[[.]], lab[[.]]))
kst <- lapply(kst, \(.) {.$clust <- factor(.$clust, exclude="x"); .})
names(kst) <- names(lab); kst$tum <- jst$tum
sapply(kst, \(.) table(.$clust, exclude=NULL))$str
BEC epi fib FRCcts FRCpv FRCtcz LEC mCAF
1034 212 653 941 744 3753 945 2016
$mye
AML DC mac.pi mono.c mono.nc pDC TAM tum
2943 1700 1529 1276 213 1587 866 1750
$nkt
NK NK/ILC Tcm Tcn Tex Tfh Tha Thn Tp Treg
42 422 197 747 1969 921 310 3513 358 1071
$tum
tum
60121
kid <- lapply(kst, \(.) .$clust)
idx <- names(kid <- unlist(unname(kid)))
kid <- kid[match(colnames(sce), idx)]
pbs <- aggregateAcrossCells(sce, kid[idx])
pbs <- logNormCounts(pbs)Code
n <- 7
ks <- colnames(pbs)
es <- logcounts(pbs)
gs <- lapply(ks, \(i) {
j <- setdiff(ks, i)
fc <- apply(es[, j], 2, \(.) es[, i]/.)
fc <- rowMins(fc)
names(tail(sort(fc), n))
}) |> unlist()
mx <- .z(t(es[gs, ]))
pal <- setNames(.pal_kid[seq_along(ks)], ks)
pheatmap::pheatmap(mx, annotation_legend=FALSE,
legend=FALSE, show_colnames=TRUE,
treeheight_row=5, treeheight_col=5,
fontsize=8, cellwidth=8, cellheight=8,
cluster_rows=FALSE, cluster_cols=FALSE,
color=pals::coolwarm(), border_color=NA,
gaps_col=cumsum(rep(n, length(ks))),
annotation_colors=list(x=pal),
annotation_names_col=FALSE,
annotation_col=data.frame(
row.names=gs,
x=rep(colnames(pbs), each=n)))appendix
saving
saveRDS(jst, "outs/jst.rds")
saveRDS(kst, "outs/lv2.rds")
session
sessionInfo()R version 4.5.1 (2025-06-13)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.6.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/Madrid
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] scran_1.36.0 patchwork_1.3.2
[3] ggplot2_4.0.0 tidyr_1.3.1
[5] dplyr_1.1.4 BiocParallel_1.42.2
[7] scuttle_1.18.0 SingleCellExperiment_1.30.1
[9] SummarizedExperiment_1.38.1 Biobase_2.68.0
[11] GenomicRanges_1.60.0 GenomeInfoDb_1.44.3
[13] IRanges_2.42.0 S4Vectors_0.46.0
[15] BiocGenerics_0.54.1 generics_0.1.4
[17] MatrixGenerics_1.20.0 matrixStats_1.5.0
loaded via a namespace (and not attached):
[1] gtable_0.3.6 xfun_0.54 htmlwidgets_1.6.4
[4] lattice_0.22-7 vctrs_0.6.5 tools_4.5.1
[7] parallel_4.5.1 tibble_3.3.0 cluster_2.1.8.1
[10] BiocNeighbors_2.2.0 pkgconfig_2.0.3 pheatmap_1.0.13
[13] Matrix_1.7-4 RColorBrewer_1.1-3 dqrng_0.4.1
[16] S7_0.2.0 pals_1.10 lifecycle_1.0.4
[19] GenomeInfoDbData_1.2.14 compiler_4.5.1 farver_2.1.2
[22] statmod_1.5.1 bluster_1.18.0 mapproj_1.2.12
[25] codetools_0.2-20 htmltools_0.5.8.1 maps_3.4.3
[28] yaml_2.3.10 pillar_1.11.1 crayon_1.5.3
[31] limma_3.64.3 DelayedArray_0.34.1 abind_1.4-8
[34] metapod_1.16.0 locfit_1.5-9.12 rsvd_1.0.5
[37] tidyselect_1.2.1 digest_0.6.37 BiocSingular_1.24.0
[40] purrr_1.1.0 labeling_0.4.3 fastmap_1.2.0
[43] grid_4.5.1 colorspace_2.1-2 cli_3.6.5
[46] SparseArray_1.8.1 magrittr_2.0.4 S4Arrays_1.8.1
[49] dichromat_2.0-0.1 edgeR_4.6.3 withr_3.0.2
[52] UCSC.utils_1.4.0 scales_1.4.0 rmarkdown_2.30
[55] XVector_0.48.0 httr_1.4.7 igraph_2.2.1
[58] ScaledMatrix_1.16.0 beachmat_2.24.0 evaluate_1.0.5
[61] knitr_1.50 irlba_2.3.5.1 rlang_1.1.6
[64] Rcpp_1.1.0 glue_1.8.0 rstudioapi_0.17.1
[67] jsonlite_2.0.0 R6_2.6.1